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Abstract 



Many aspects of the historical relationships between populations in a species are reflected in genetic 
data. Inferring these relationships from genetic data, however, remains a challenging task. In this 
paper, we present a statistical model for inferring the patterns of population splits and mixtures 
in multiple populations. In this model, the sampled populations in a species are related to their 
common ancestor through a graph of ancestral populations. Using genome-wide allele frequency 
data and a Gaussian approximation to genetic drift, we infer the structure of this graph. We 
applied this method to a set of 55 human populations and a set of 82 dog breeds and wild canids. 
In both species, we show that a simple bifurcating tree does not fully describe the data; in contrast, 
we infer many migration events. While some of the migration events that we find have been 
detected previously, many have not. For example, in the human data we infer that Cambodians 
trace approximately 16% of their ancestry to a population ancestral to other extant East Asian 
populations. In the dog data, we infer that both the boxer and basenji trace a considerable fraction 
of their ancestry (9% and 25%, respectively) to wolves subsequent to domestication, and that 
East Asian toy breeds (the Shih Tzu and the Pekingese) result from admixture between modern 
toy breeds and "ancient" Asian breeds. Software implementing the model described here, called 
TreeMix, is available at http://treemix.googlecode.com. 



Author Summary 



With modern genotyping technology, it is now possible to obtain large amounts of genetic data 
from many populations in a species. An important question that can be addressed with these 
data is: what is the history of these populations? There is a long history in population genetics 
of inferring the relationships among population as a bifurcating tree, analogous to phylogenetic 
trees for representing the evolution of species. However, it has long been recognized that, since 
populations from the same species exchange genes, simple bifurcating trees may be an incorrect 
representation of population histories. We have developed a method to address this issue, using a 
model which allows for both population splits and gene flow. In application to humans, we show 
that we are able to identify a number of both previously known and unknown episodes of gene 
flow in history, including gene flow into Cambodia of a population only distantly related to modern 
East Asia. In application to dogs, we show that the boxer and basenji breeds have a considerable 
component of ancestry from grey wolves subsequent to domestication. 

Introduction 

The extant populations in a species result from an often-complex demographic history, involving 
population splits, gene flow, and changes in population size. It has long been recognized that 
genetic data can be used to learn about this history [1-3]. In humans, early approaches to inferring 
history from genetics were limited to using a relatively small number of blood group or other protein 
polymorphisms [1, 4-6]. These types of studies were then superseded by analyses of DNA markers, 
which have progressed from single marker studies [3] to studies involving hundreds of thousands of 
markers [7]. It is now feasible to collect genome-wide genetic data in any species; to a large extent 
it is no longer the data collection, but rather the statistical models used for analysis, that limit the 
historical insight possible. 

There are many statistical approaches to demographic inference from genetic data. One ap- 
proach is to develop an explicit population genetic model for the history of a set of populations, 
framed in terms of the effective population sizes of the populations, the times of population splits, 
the times of demographic events (such as population bottlenecks), and other relevant parameters. 
The values of these parameters can then be learned from the data using a variety of techniques, 
often involving simulation [8-16]. These approaches have the advantage of allowing flexible model- 
ing of a wide variety of demographic scenarios, but the disadvantage that they can only be applied 
to one or a few populations at a time. 

Another type of approach to learning about population history uses methods that summarize the 
major components of genetic variation in a sample by clustering or principal components analysis 
[17-20]. Although these methods do not model history explicitly, the inferred components can 
often be interpreted post hoc as representing historical populations, and individuals or populations 
that are mixtures of different components as evidence of admixture between these populations 
(e.g., [17, 21-23]). However, these methods are not directly informative about history; indeed, the 
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relationship between the major components of genetic variation and true underlying demography 
is not always intuitive [24-26]. 

A different class of approaches focuses on the relationships between populations, by representing 
a set of populations as a bifurcating tree [1, 27-32]. In these models, the details of the demographic 
histories of the population are absorbed into the branch lengths of the tree [1, 33]. This approach 
has the advantage of being applicable to large numbers of populations; however, a major caveat 
when modeling the history of populations as a tree is that gene flow violates the assumptions of the 
model [2, 34, 35]. It is often difficult to know, a priori, how well the history fits a simple bifurcating 
tree. Explicit tests for the violation of a tree model have been developed [35-40]. These tests have 
been used, most notably, to infer the existence of gene flow between modern and archaic humans 
[39, 41, 42], as well as between diverged modern human populations [37, 43, 44]. 

In this paper, we present a unified statistical framework for building population trees and testing 
for the presence of gene flow between diverged populations. In this framework, the relationship 
between populations is represented as a graph, allowing us to model both population splits and 
gene flow. Graph-based models are of growing interest in phylogenetics [45, 46], but have been 
rarely used in population genetics (with some exceptions [37, 40, 47]). 

Results 

The starting point for our model was first proposed by Cavalli-Sforza and Edwards [1], and we draw 
additionally on related models by Nicholson et al. [33] and Coop et al. [48]. Our goal is to provide a 
statistical framework for inferring population networks that is motivated by an explicit population 
genetic model, but sufficiently abstract to be computationally feasible for genome-wide data from 
many populations (say, 10-100). Our primary aim is to represent the topology of relationships 
between populations, rather than the precise times of demographic events. 

Our approach to this problem is to first build a maximum likelihood tree of populations. We 
then identify populations that are poor fits to the tree model, and model migration events involving 
these populations. Below, we first describe this approach in an idealized setting, and then describe 
the modifications necessary for implementation in practice. 

Model 

In the most simple case, consider a single SNP, and let the allele frequency of one of the alleles at 
this SNP in an ancestral population be xa (We use a lowercase x to denote that this is a parameter 
rather than a random variable. We initially consider distributions conditional otixa)- Now consider 
a descendant population B. We model X B , the allele frequency of the SNP in population B, as: 

X B = x A + ee (1) 

with 

e B ~ N{0,c B x A [l-x A ]) (2) 
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A. Example tree B. Covariance matrix for tree in A. 




1 1 1 1 1 1 X, X 2 X 3 X 4 

0.00 0.05 0.10 0.15 0.20 0.25 



Drift parameter 



C. Example graph D. Covariance matrix for graph in C. 




' 1 1 1 1 1 X, X 2 X 3 X 4 

0.00 0.05 0.10 0.15 0.20 0.25 



Drift parameter 

Figure 1: Simple examples. A. An example tree. B. The covariance matrix implied by the tree 
structure in A. Note that the covariance here is with respect to the allele frequency at the root, 
and that each entry has been divided by xa[1 — xa] to simplify the presentation. C. An example 
graph. The migration edge is colored red. Parental populations for population 3 are labeled Y 
and Z; see the main text for details. D. The covariance matrix implied by the graph in C; again, 
each entry has been divided by xa[1 — xa]- The migration terms are in red, and the non-migration 
terms are in blue. 
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where eg is a factor that corresponds to the amount of genetic drift that has occurred between 
the ancestral population and B. This Gaussian model was first introduced by Cavalli-Sforza and 
Edwards [1], and the motivation for this model is outlined in Nicholson et al. [33]. Briefly, if the 
amount of genetic drift between the two populations is small (at most on a timescale of the same 
order as the effective population size), then the diffusion approximation to a Wright-Fisher model 
of genetic drift leads to Equation 2 with c B ~ where t is the number of generations separating 
the two populations, and N e is the effective population size [33]. We do not model the boundaries 
of the allele frequencies at zero and one, nor do we consider new mutations. This means that 
this model will be most accurate for alleles that were at intermediate frequency in the ancestral 
population. 

Now consider a descendant population of B; let us call this population C, and the allele fre- 
quency in the population Xc- Using the same model: 

X c = X B + e c (3) 
= x A + e B + e c (4) 

where 

ec~N(0,c c X B [l-X B }). (5) 
We can write down the expectation and variance of Xc as: 

E[X C \ = E[x A + e B + e c ] (6) 

= x A (7) 

and: 

Var(Xc) = Var(x A + e B + e c ) (8) 

= Var{e B ) + Var(e c ) + 2Cov{e B , e c ). (9) 

We then assume that the amount of genetic drift between all the populations is small. This implies 
that X B [l — X B ] is well-approximated by xa[1 — xa\ in Equation 5, and hence the amount of genetic 
drift between A and B is approximately independent of the amount of genetic drift between B and 
C [35]. With these simplifications: 

Var(Xc) ~ Var(e B ) + Var(e c ) (10) 
~ (c B + c c )x A [l - x A ]- (11) 

We thus have a model for Xc, conditional on xa- 

X c ~N(x a ,(cb + cc)xa[1-x a ]). (12) 
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Multiple populations. Now consider a set of four populations, all related to an ancestral 
population by a tree, as depicted in Figure 1A. Let the allele frequencies in the four popula- 
tions be denoted Xi, X2, A3, and A4, respectively, and the vector of all four frequencies be X. 
We want to write down a joint distribution for X given the tree. We start by writing down 
the covariance between any two populations with respect to the ancestral allele frequency (i.e. 
Cov(Xi, Xj) = E[(Xi — xa){Xj — xa)})- This is simply the variance of the common ancestor of the 
two populations: 

Cov(X 1 ,X 2 ) = c 2 x a [1-xa] (13) 
Cov{X^X i ) = c l x A [l-x A ] (14) 
Cov(X 1 ,X 3 ) = (15) 

and so on (Figure IB). 

Let us use V to denote the variance-covariance matrix of allele frequencies between populations 
implied by the tree. Now, if we knew the value of xa, we could model X as: 

X ~ MVN(xa, V) (16) 

where xa = [xa,xa,%a,xa] and MVN denotes the multivariate normal distribution. The covari- 
ance matrix with respect to the root, V, is distinct from the sample covariance matrix; we return 
to this distinction later. This multivariate normal model was proposed by Felsenstein [28]. Addi- 
tionally, a multivariate normal model was used by Coop et al. [48] and Weir and Hill [49], although 
these authors did not explicitly model the variance-covariance matrix in terms of a tree. 

Modeling migration. To extend this framework to include migration, we allow populations 
to have ancestry from multiple parental populations [35, 40]. The contribution of each parental 
population is weighted; if we assume admixture occurs in a single generation, these weights are 
related to (though not necessarily equivalent to) the fraction of alleles in the descendant population 
that originated in each parental population. Consider population 3 in Figure 1C (recall that the 
allele frequency in this population is A3). We have labeled the two parental populations Y and Z; 
let the allele frequencies in these populations be Ay and Xz, respectively. We model A3 as: 

X 3 = wX Y + (l-w)(X z + e 3 ) (17) 

where €3 ~ N(0, C3Xa[1 — xa])- Note that there is some question as to how to weight €3, the genetic 
drift specific to population 3. In reality, it comes from three sources: drift since Y but before the 
population mixture, drift since Z but before the population mixture, and drift since the mixture. 
These three components should have weights w, 1 — w, and 1, respectively. However, the three 
components are not all separately identifiable. For ease of computation, we estimate only a single 
drift parameter in this situation, and weight it by 1 — w (Supplementary Information). 

Additionally, there is a choice of whether the edge from Z or the edge from Y should be 
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considered the "migration" edge; these two possibilities not identifiable in the model. In practice, 
we set the edge with the largest weight as the non- migration edge, and the other edge (or edges) 
as the "migration" edge(s). 

With these simplifications, the variance of A3 can be written in the mixture case as: 

Var(X 3 ) = Var(wX Y + (1 - w)(X z + e 3 )) (18) 
= w 2 Var(X Y ) + (1 - w) 2 [Var{X z ) + Var(e 3 )} + 2w(l - w)Cov(X Y ,X z ) (19) 

We can now consider multiple populations related by a graph instead of a tree (Figure 1C). The 
variance-covariance matrix V can be filled in as before, but now including terms for migration 
(Figure ID). This model can be written in terms of a directed acyclic graph (the lack of cycles 
follows from the fact that no population can contribute genetic material to its own ancestor), 
where the c parameters correspond to edge lengths (Supplementary Material). For subsets of up 
to four populations, this model is closely related to the "/— statistics" used as tests for treeness 
by Reich et al. [37] (Supplementary Material). 

Normalization. As described above, V depends on the ancestral allele frequency xa- This means 
that the true variance-covariance matrix will differ by a scaling factor between SNPs with different 
values of xa- In much work on this type of model, investigators have normalized allele frequencies 
to account for this. One potential normalization is the arcsine square-root transformation [1]. 
However, a drawback to this normalization is that it is non-linear; the expected value of the allele 
frequency in the descendant populations is no longer xa, but pushed towards the boundaries, which 
could induce spurious correlations between the most drifted populations [50]. Another plausible 
transformation would be to scale all allele frequencies by — ft), where /t is the mean allele 
frequency across populations [19, 33]. Both of these transformations increase the influence of 
polymorphisms that were rare in the ancestral population. However, these are precisely the loci 
where the approximation of our model to a true population genetics model is most likely to break 
down. For these reasons, we choose to work directly with untransformed allele frequencies. 

Properties of the sample covariance. In practice, the multivariate normal model in Equation 
16 is impractical because we do not know the ancestral values of allele frequencies, but instead only 
the values in sampled descendant populations. This means that V cannot be estimated directly from 
data. However, consider instead the sample covariance matrix W, where Wjj = E[{Xi— jl){Xj— £1)], 
where ft = ^ — l ' m 1S ^ ne numDer °f populations, and Aj and Xj are the sample allele frequencies 
in populations i and j. W is closely related to V as follows: 
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Wij = E[(Xi - ^(Xj - jl)} (20) 
= E[(Xi - x A - A + x A )(Xj - x A - A + (21) 
= - x A )(Xj - x A ) - (Xi - x A )(fi - x A ) ~ (Xj - x A )(fi - x A ) + (A - x A f] (22) 

1 m .. m 1 m m 

= v ^ - ™ E v ^ - - E v ^ + -2 E E ( 23 ) 
^=1 fc=i fc=ifc'=i 

In the following section, we will describe how we perform inference based on the sample covariance 
matrix W. 

Finite samples and multiple (potentially correlated) SNPs. Now assume that we have 
genotyped n SNPs in each of m populations. Let the sample allele frequency at SNP k in population 
i be Xik- We can estimate the sample covariance matrix W: 

W = *^2=l\( X ik ~ fik)(X jk ~ Afc)] / 24 n 

13 ' n 

where Afc = ^ YaLi ^ik- The fact that in practice we have finite samples from a population has 
two effects on this matrix. First, sampling leads to a biased estimation of the terms on the diagonal, 
since sampling can be thought of as adding an amount of "drift" to each population (as well as 
smaller effect on the off-diagonal terms; see Supplementary Material for details). Additionally, it 
leads to some uncertainty in the off-diagonal terms. To account for the biased diagonal terms, in 
practice we calculate a corrected version of W that removes this bias (Supplementary Material). 
To account for uncertainty in the values of this matrix, we use a block resampling approach (see 
below). Finally, with multiple SNPs, we are working with SNPs with many different values of xa- 
In this case, the x a[1 — xa] terms described above can be thought of as xa\\ — xa]', he., the mean 
across SNPs of xa\\ — xa]- 

We now want to write down a likelihood for W given W. One possibility would be to use the 
Wishart distribution, since the sample covariance matrix of multivariate normal random variables 
has this form. However, computation of the Wishart density involves computationally-intensive 
matrix inversion, so we took an alternative approach. Consider the observed covariance between 
populations i and j, Wjj. If we had a large number of independent genomic regions and estimated 
this quantity separately in each independent region, the sampling distribution would be approx- 
imately normal with mean Wjj (by appeal to the central limit theorem). We thus model Wjj 
as: 

W -~JV(Wtf,<7?-) (25) 

where Oij is the standard error in the estimation of Wy. Because the allele frequencies at nearby 
SNPs are correlated (i.e., there is linkage disequilibrium), a naive estimate of <jjj that treated each 
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SNP as independent would be too small. We instead take a resampling approach to estimate aij. 
We split the genome into p blocks, such that there are K SNPs per block (with K chosen so that 
the block sizes are larger than blocks of linkage disequilibrium) [36]. (If K does not divide evenly 
into n, we discard the remaining SNPs.) We then calculate W separately in each block. Let Wyfc 
be the sample covariance between two populations i and j in block k. Now, 

W (J = ^ k=l ™ l]k (26) 

and 



If we take each pair of populations in turn, we can write down a composite likelihood for W: 

m m 

l(wiw) -nn^w^-^) ( 2§ ) 
i=i j=i 

where N(Wij\Wij, cr|-) is a Gaussian density with mean Wjj and variance afj evaluated at Wjj. 

Finally, we wanted to define measures for how well the model fits the data. First, we define the 
matrix of residuals in this model, R. These quantities are useful for visualization and fitting: 

R = W-W. (29) 

Positive residuals indicate pairs of populations where the model underestimates the observed covari- 
ance, and thus populations where the fit might be improved by adding additional edges. Negative 
residuals indicate pairs of populations where the model overestimates the observed covariance; 
these are a necessary outcome of having positive residuals, but can also sometimes be interpreted 
as populations that are forced too close together due to unmodeled migration elsewhere in the 
graph. These residuals can be used to define a measure of the fraction of the variance in W that 
is explained by W. Let us call this fraction /: 

Em v^m i) _ tj-n.2 

, ! i=l 2^j=i+l K^ij *V , 3Q s 

££i£r=i + i(w y 

where R = i= m{m-t)/2 ~ and ^ = i= m(J-i)/2 13 ■ T ^ is q uan tity approximates the fraction of 
the variance in relatedness between populations that is accounted for by the model. 



Estimation. We implemented an algorithm, called TreeMix, that searches for the graph that 
maximizes the composite likelihood in Equation 28. A search that enumerates all graphs is infeasible 
unless m is very small, so to simplify the search we make the assumption that the history of the 
sampled populations is approximately tree-like. We thus start by searching for the maximum 
likelihood tree, taking an algorithmic approach similar to Felsenstein [30]. 
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After building the tree, we fix the position of the root. (In the tree model the position of the 
root is not identifiable, as the evolution of allele frequencies along the tree is reversible under the 
Gaussian model when drift is assumed to be small. In a graph model, though the position of the 
root is partially identifiable, in all applications we assume that the position of the root is fixed using 
prior information about known outgroups). We then calculate the residual covariance matrix, R, 
and add migration edges in a directed matter. First, we find the M pairs of populations with the 
maximum residuals. We then attempt adding a migration edge between populations in the vicinity 
of each of the M population pairs. For each attempted graph (or tree) topology, we optimize the 
branch lengths and migration edge weights (Methods). 

After finding the single migration edge that most increases the likelihood, we attempt a series 
of local changes to the graph structure (Methods). We then iterate over this procedure to add 
additional migration edges. In principle, migration edges could be added until they are no longer 
statistically significant (see the following paragraph). In our experience in practice, however, we 
prefer to stop adding migration events well before this point so that the resulting graph remains 
interpretable. 

Significance testing. After building the maximum likelihood graph, we would like to quantify 
our uncertainty in the resulting graph structure. In particular, we would like to quantify our 
confidence in individual migration events. However, because the likelihood in Equation 28 is a 
composite likelihood, it cannot be used directly for formal tests for significance. Instead, we take 
a resampling approach to test the support for individual migration edges. 

Consider a given migration edge, with corresponding weight w. We wish to calculate a p-value 
for this weight (under the null hypothesis that w = 0, and for a fixed graph structure). To do this, 
we use the Wald statistic — rW, where se(w) is the standard error in the estimate of the weight, 

se(w) ' v / o > 

which is distributed iV(0, 1) under the null. To obtain the standard error, recall that we have split 
the genome into p independent blocks. We use the jackknife estimates of both w and the standard 
error in w (where we jackknife over blocks). Let i index blocks, and w.{ be the estimated weight 
computed using all blocks except i. Then: 




This allows us to calculate a p-value for the migration edge. There are a number of complications 
to the interpretation of this p-value. First, there is the issue of multiple testing-there are at least 
2m— 2 edges in the graph (recall that m is the number of populations), and thus around 4m 2 possible 
migration events. More importantly, the p-value is generated under a heavily parameterized model: 
we are comparing a fixed graph structure with a migration event to that same graph without the 
migration event. A "significant" p-value simply indicates that the hypothesized migration event 
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significantly improves the fit to the data; this does not account for the possibility of errors in the 
graph structure, or indicate that the particular migration event tested is the correct one (rather 
than a migration event between a different pair of populations). For this reason, we treat the 
precise p-value generated by this procedure with caution, and use additional, less-parameterized 
methods like three- and four-population tests [37] to test the robustness of the inference. 

Simulations 

We tested the performance of the TreeMix method in simulations. We generated coalescent simu- 
lations from several histories; the basic structure was a set of 20 populations produced by a serial 
bottleneck model like that used by DeGiorgio et al. [51] to model human history (Figure 2A). The 
parameters of the simulations were chosen to be reasonable for non- African human populations; 
we used an effective population size of 10,000, and a history where all 20 populations share a 
common ancestor 2000 generations in the past. Each individual simulation involved 400 regions of 
approximately 500kb each, and thus recapitulated many aspects of real data, including hundreds 
of thousands of loci and the presence of linkage disequilibrium. 

Tree simulations. First, we tested the performance of the algorithm on truly tree-like data. 
We generated 100 independent simulations of 20 chromosomes from each population using the 
above demographic model without migration, and inferred population trees. The inferred trees 
perfectly matched the simulated model in all cases (Figure 2B, Supplementary Figure 2), and the 
fitted tree model accounted for an average of 99.8% of the variance in W. To test the effect of 
SNP ascertainment, we then inferred trees using only SNPs that were polymorphic in one of the 
populations (either population 1 or population 20); this ascertainment scheme did not decrease 
accuracy of the inferred topology, though it did alter the inferred branch lengths (Supplementary 
Figure 3). 

We used these simulations without migration to test the calibration of our p- values for migration 
events. For each simulation, after building the maximum likelihood tree, we introduced a migration 
event between two random populations and tested it for significance. As expected if the p-values 
are properly calibrated, their distribution is approximately uniform (Supplementary Figure 4). 

Finally, we performed tree simulations in a situation where fixed differences and new mutations 
(rather than shared polymorphisms inherited from a common ancestor) were common between the 
populations; in this context the population genetic interpretation of the model breaks down. We 
performed simulations where all the true branch lengths were 50 times longer than in the original 
model, corresponding to a history where the 20 populations share a common ancestor approximately 
100,000 generations in the past. Again, we see no errors in the topology of the inferred trees 
(Supplementary Figures 5, 6). In this situation, the covariances between closely-related populations 
tend to be slightly underestimated; in more extreme situations this could lead to spurious inferences 
of migration (Supplementary Figures 5, 6). However, overall, these simulations suggest that the 
model will still be useful even in situations where the population genetic interpretation is not strictly 
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correct. 



Simulations with migration. We then introduced migration events into our simulations. We 
generated simulations under the same model described above; however, we now simulated an ad- 
mixture event approximately 100 generations before the present where one population receives a 
fraction of its ancestry (either 10% or 30%) from one of the other populations. We tried ten dif- 
ferent pairwise combinations of populations, and generated 100 simulations for each pair. For each 
simulation, we ran TreeMix and allowed it to infer a migration event. We then judged the error 
rate of the algorithm as the fraction of times the inferred topology of the graph was not exactly 
correct (this is a conservative estimate of the error rate, in that inferred graph topologies that are 
very close to the truth are counted as errors). In general, TreeMix was able to correctly infer the 
graph structure in these simulations (Figure 2C). However, accuracy dropped considerably when 
migration was between closely related populations without outgroups present in the data (these 
are populations 1 and 20 in the model; Figure 2C). The major types of errors produced in the 
simulations were incorrectly inferred directions of migration arrows and inference of admixture 
in populations related to the truly admixed population (Supplementary Material, Supplementary 
Figure 7). 

We next asked whether the mixture "weights" inferred in the model can be interpreted as 
admixture proportions. To do this, we simulated admixture events of varying proportion between 
the first and tenth population in the serial bottleneck model described above, set the graph to the 
true topology, and estimated the mixture weight. The weights are indeed correlated with the true 
ancestry fraction, but underestimate relatively high admixture proportions in these simulations 
(Figure 2D). The precise bias in the estimation of this parameter will depend, in real data, on 
largely unknowable parameters (Supplementary Material). 

Application to humans 

To test the performance of the TreeMix model with real data, we applied it to humans, whose genetic 
history has been studied extensively [7, 21, 52, 53]. We applied the model to a dataset consisting of 
about 125,000 SNPs ascertained by low-coverage genome sequencing in a single Yoruban individual 
and then genotyped in 55 modern and archaic human populations [54]. In all that follows, we 
excluded the two Oceanian populations because they gave inconsistent results across datasets. We 
believe this difficulty results from the fact that these populations contain ancestry from multiple 
sources, making the graph estimation somewhat unstable when they are included (Supplementary 
Material). We first built the tree of all 53 remaining populations (Figure 3A). This tree largely 
recapitulates the known relationships among population groups [7], and explains 98.8% of the 
variance in relatedness between populations (though this is high, it is less than the 99.8% observed 
in the simulations of a true tree model). We examined the residuals of the model's fit to identify 
aspects of ancestry not captured by the tree (Figure 3B). A number of known admixed populations 
stand out: in particular, these include the Mozabite and Middle Eastern populations. 



11 



A. Simulated model B. Inferred tree 




0.000 0.004 0.008 



Drift parameter 



C. Accuracy of graph inference D. Admixture weight estimation 




source 1 1 1 1 5 5 5 10 10 15 1 1 1 1 1 1 

dest. 51015 20 1015 20 15 20 20 qq q-| q 2 03 04 05 

Simulated migration event 

Simulated weight 



Figure 2: Performance on simulated data. A. The basic outline of the demographic model 
used. B. Trees inferred by TreeMix. We simulated 100 independent data sets, under the 
demographic model in A., and inferred the tree. All simulations gave the same topology; plotted 
are the mean branch lengths. C. Performance in the presence of migration. We added 
migration events to the tree in A. and inferred the structure of the graph. Each point represents 
the error rate over 100 independent simulations, defined as the fraction of simulations where the 
inferred graph topology does not perfectly match the simulated topology. On the x-axis we show 
the populations involved in the simulated migration event; e.g., if the source population is 1 and 
the destination population is 10, this is a migration event from population 1 to population 10, as 
labeled in A. D. Admixture weight estimation. We simulated admixture events with different 
weights from population 1 to population 10, and inferred the weight. Each point is the mean across 
100 simulations, and the bar represents the rang^ 



We then sequentially added migration events to the tree. In Figure 4, we show the inferred 
graph with ten migration edges; p-values for all reported migration edges are less than 1 x 1CT 30 
(we show the graph with the maximum likelihood over several independent runs of TreeMix with 
random orders of input populations). This graph model explains 99.8% of the variance in relat- 
edness between populations. As expected from examination of Figure 3B, the migration events 
recapitulate many known events in human history. We infer that the Mozabite are the result of 
admixture between an African and a Middle Eastern population (with about 33% of their ancestry 
from Africa), and that Middle Eastern populations also have African ancestry (Palestinians and 
Bedouins: w = 13% from Africa; Druze: w = 6%). This is consistent with previously reported ad- 
mixture proportions from these populations [43, 55]. Additionally, we identify the known European 
ancestry in the Maya (w = 12%) [21], and infer that the Uyghur and Hazara populations are the 
result of admixture between west Eurasian and East Asian populations (w = 46% and 47% from 
west Eurasia, respectively) [20, 21, 56]. 

Several additional migration events in the human data have not been previously examined in 
detail, but are consistent with previous clustering analysis of these populations [7, 20, 21]. These 
include migration from Africa to the Makrani and Brahui in Central Asia (w = 5%) and from a 
population related to East Asians and Native Americans (which we interpret as likely Siberian) to 
Russia (w = 11%). 

Two inferred edges were unexpected. First, perhaps the most surprising inference is that Cam- 
bodians trace about 16% of their ancestry to a population equally related to both Europeans and 
other East Asians (while the remaining 84% of their ancestry is related to other southeast Asians). 
This is partially consistent with clustering analyses, which indicate shared ancestry between Cam- 
bodians and central Asian populations [7]. To confirm that the Cambodians are admixed, we 
turned to less parameterized models. The predicted admixture event implies that allele frequencies 
in Cambodia are more similar to those in African populations than would be expected based on 
their East Asian ancestry. To test this, we used three-population tests [37]. We tested the trees 
[African, [Cambodian, Dai]] for evidence of admixture in the Cambodians (Methods). When using 
any African population, there is strong evidence of admixture (when using Yoruba, Z = —7.0 
[p = I x 10~ 12 ]; when using Mandenka, Z = -7.3 [p = 1 x 10~ 12 ]; when using San, Z = -4.8 
[p = 8x 10~ 7 ]). We conclude that the Cambodian population is the result of an admixture event 
involving a southeast Asian population related to the Dai and a Eurasian population only distantly 
related to those present in these data. 

Finally, we infer an admixture edge from the Middle East (a population related to the Mozabite, 
a Berber population from northern Africa) to southern European populations (w = 22%). This 
migration edge is the one edge that is not consistent across independent runs of TreeMix on these 
data (Supplementary Figure 8). In particular, an alternative graph (albeit with lower likelihood) 
places the Mozabite as an admixture between southern Europe and Africa (rather than the Middle 
East and Africa), and does not include an edge from the middle East to southern Europe. We thus 
hesitate to interpret this result, except to note that the relationship between northern African, 
the Middle East, and southern Europe involves complex patterns of gene flow that merit further 
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investigation [43, 57]. 

To test the robustness of our results to SNP ascertainment, we additionally ran TreeMix on the 
same set of populations using a set of SNPs ascertained in a single French individual. The inferred 
graph was nearly identical (Supplementary Figure 10). Additionally, as noted above, different 
random input orders for the populations gave very similar results (Supplementary Figure 8). We 
conclude from this that the model is able to consistently and accurately infer the major mixture 
events in the history of a species. This approach is computationally efficient: building the tree took 
around five minutes on a standard desktop computer (with a processor speed of 3.1 GHz), and 
adding ten migration events to the tree took about four and a half hours (the major computational 
cost is in the iterative estimation of migration weights) . 

Application to dogs 

While human populations have been extensively studied, we next applied the model to dogs, a 
species where considerably less is known about population history. In particular, we applied the 
model to a dataset consisting of about 60,000 SNPs genotyped in 82 dog breeds or wild canids 
[58]. As for humans, we first inferred the maximum likelihood tree (Figure 5A). The differences in 
history between dogs and humans are striking: there are long terminal branches leading to each dog 
breed in the inferred tree (Figure 5A, recall that the terminal branch lengths account for sample 
size). This is consistent with the known strong bottlenecks in the establishment of dog breeds 
[23]. However, examining the residuals from the model revealed a number of populations that do 
not fit a strict tree model (Figure 5B); indeed, the tree model explained 94.7% of the variance in 
relatedness between breeds, somewhat less than between human populations. 

We sequentially added migration events to the tree in Figure 5A. In Figure 6, we show the in- 
ferred graph with ten migration events, which explains 96.8% of the variance in relatedness between 
breeds (which suggests that additional events exist in the data). In the following paragraphs, we 
describe some of these events. 

We infer that the bull mastiff is the result of an admixture event between bulldogs and mastiffs. 
This is a known event [59] ; we estimate the admixture proportions as 33% bulldog and 67% mastiff. 
We further examined this event using four-population tests for treeness. As expected given the 
known history, the tree [[boxer, bulldog], [mastiff,bull mastiff]] fails the four-population test (Z = 3.5, 
p = 0.002), while replacing the bull mastiff with other related breeds that we do not predict to 
be involved in the admixture event results in trees that pass this test. For example, the tree 
[[boxer, bulldog], [mastiff,Boston terrier]] passes the four-population test with Z = —0.3. 

The most visually apparent residuals in Figure 5B are accounted for in the graph by an admix- 
ture event from the grey wolf into the basenji, an ancient African breed of dog (w = 25%). Such a 
high mixture fraction is consistent with previous clustering analyses of these data [23, 60]. We again 
sought to confirm this signal in a less-parameterized model. We tested the four-population tree 
[[wolf, ancient breed], [basenji, Afghan hound]] with various "ancient" dog breeds. We could not find 
a tree that passed the four-population test (with Akita as the ancient breed, Z = 11.7, p < 1 x 10 -30 ; 
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A. Maximum likelihood human tree 




Figure 3: Inferred human tree. A. Maximum likelihood tree. Plotted is the maximum- 
likelihood tree. Populations are colored according to geographic location (black: archaic humans, 
red: Africa, brown: Middle East, green: Europe, blue: Central Asia, purple: America, orange: 
East Asia). The scale bar shows ten times the average standard error of the entries in the sample 
covariance matrix (W). For analysis including Oceania, see Supplementary Figures 11 and 12. B. 
Residual fit. Plotted is the residual fit from the maximum likelihood tree in A. We divided the 
residual distance between each pair of populations i and j by the average standard error across all 
pairs. We then plot in each cell this scaled residual. Colors are described in the palette on 
the right. Residuals above zero represent populations that are more closely related to each other 
in the data than in the best-fit tree, and thus are candidates for admixture events. 
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Figure 4: Inferred human tree with mixture events. Plotted is the structure of the graph 
inferred by TreeMix for human populations, allowing ten migration events. Migration arrows are 
colored according to their weight. Horizontal branch lengths are proportional to the amount of 
genetic drift that has occurred on the branch. The scale bar shows ten times the average standard 
error of the entries in the sample covariance matrix (W). The residual fit from this graph is 
shown in Supplementary Figure 9. Admixture from Neandertals to non- African populations is only 
apparent when considering subsets of the data (see Discussion and Supplementary Figure 15). 
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with Alaskan Malamute, Z = 13.0, p < 1 x 1CT 30 ), confirming the presence of gene flow in these trees. 
Replacing the basenji with the saluki in these analyses resulted in trees that pass the four-population 
test (for example, the tree [[wolf, Akita] , [Afghan hound, saluki]] passes with Z = — 0.03, p = 0.51). 
Though we cannot have complete confidence in the precise migration events, these results are 
consistent with admixture between gray wolves and the basenji. 

Another breed that stands out in this analysis is the boxer (Note that many of the SNPs 
used in this study were ascertained using a boxer individual, so we may have increased power to 
identify migration events involving this breed). We infer a significant genetic contribution from 
wolves to the boxer (w = 8%), and migration between the boxer and the Chinese shar-pei, a 
distantly-related ancient breed (w = 8%). To further examine these events, we again turned to 
four-population tests. To evaluate the wolf mixture, we tested the tree [[wolf, ancient breed], [boxer, 
bulldog]]. We did not find a tree that passed the four-population test (with Akita as the ancient 
breed, Z = 3.1, p = 0.001; with Afghan Hound, Z = 3A,p = 0.0003). Replacing the Boxer with 
the Mastiff in these analyses led to trees that passed the four-population test (for example, with 
Akita as the ancient breed, Z = 0.3, p = 0.38). To evaluate the gene flow from the Boxer to 
the Chinese shar-pei, we tested the tree [[Chinese shar-pei, Akita], [boxer, bulldog]]; this tree fails 
the four-population test [Z = 3.0, p = 0.001), while the tree [[Chow Chow, Akita], [boxer, bulldog]] 
passes (Z = -0.48, p = 0.3). 

Previous analyses of these data have noted that the "toy breeds" of dog cluster together [23] . We 
find that the Chinese toy breeds (the Pekingese and the Shi Tzu) result from admixture between a 
population related to ancient East Asian dog breeds and a modern population related to the Brussels 
griffon and the pug (w = 28% from the East Asian breeds). To confirm the presence of gene flow, 
we tested four-population trees of the form [[Asian toy breed, Akita/Chow Chow], [Pug, mastiff]]. 
These trees fail, with varying levels of significance, ranging from [[Chow Chow, Shi Tzu], [Pug, 
mastiff]] (Z = -2.7, p = 0.003) to [[Akita, Pekingese], [Pug, mastiff]] (Z = -4.7, p = 1 x 10~ 6 ). 

Finally, we noticed that two of the sighthounds (the Borzoi and the Italian greyhound) do not 
cluster with the other sight hounds in the tree, namely greyhound, whippet and Irish wolfhound 
(Figure 5A); however, they do show evidence of having sighthound admixture in the graph (Figure 
6). These are the borzoi (which appear to be admixed between an ancient or spitz-breed dog, with 
47% ancestry from the sighthounds) and the Italian Greyhound (which appears to be admixed 
with a toy breed, with 34% ancestry from the sighthounds). This is consistent with the known 
phenotypic characteristics of these dogs; the borzoi is considered a saluki-like breed, and the Italian 
greyhound is phenotypically a small version of a greyhound [59]. 

Overall, we conclude that there has been considerable gene flow between dog breeds over the 
course of domestication; there are many additional migration events that merit further examination 
(Figure 6, Supplementary Material). 
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A. Maximum likelihood dog tree 
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B. Residual fit from tree 




Figure 5: Inferred dog tree. A. Maximum likelihood tree. Populations are colored according 
to breed type. Dark blue: wild canids, grey: ancient breeds, brown: spitz breeds, black: toy dogs, 
red: spaniels, maroon: scent hounds, dark red: working dogs, light green: herding dogs, light blue: 
mastiff-like dogs, purple: small terriers, orange: retrievers, dark green: sight hounds. The scale bar 
shows ten times the average standard error of the entries in the sample covariance matrix (W). B. 
Residual fit. Plotted is the residual fit from the maximum likelihood tree in A. We divided the 
residual distance between each pair of populations i and j by the average standard error across all 
pairs. We then plot in each cell [i,j] this scaled residual. Colors are described in the palette on the 
right. 
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Drift parameter 



Figure 6: Inferred dog graph. Plotted is the structure of the graph inferred by TreeMix for dog 
populations, allowing ten migration events. Migration arrows are colored according to their weight. 
The scale bar shows ten times the average standard error of the entries in the sample covariance 
matrix (W). See the main text for discussion. The residual fit from this graph is presented in 
Supplementary Figure 13. 
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Discussion 



In this paper, we have developed a unified model for inferring patterns of population splits and 
mixtures from genome-wide allele frequency data. We have shown that this model is accurate in 
simulations, largely recapitulates the known relationships between well-studied human populations, 
and is able to identify new relationships between populations in both humans and dogs. 

The TreeMix model can be thought of as a complement to methods for the identification of 
population structure [18-20] . These latter methods are powerful tools for clustering together indi- 
viduals into relatively homogenous populations (and to identify individuals that are genetic outliers 
in their population) [18-20]. However, once population structure in a species has been identified, 
these methods are not well-suited for describing how it arose, and are only indirectly informative 
about the historical relationships between different populations. The model developed in this paper 
is designed to more directly address these historical questions. 

Modeling assumptions. There are a number of assumptions, both implicit and explicit, in the 
interpretation of the TreeMix model. First, we have motivated the model in terms of inferring 
the historical splits and mixtures of populations. However, a given covariance structure of allele 
frequencies between populations can be a consequence of either a non-equilibrium demography 
(population splits and mixtures) or an equilibrium demography (populations at long-term stasis 
with a fixed migration structure) [2]. For the species analyzed in this paper, population equilibrium 
over the entire species range is not a tenable hypothesis; however, some subsets of populations may 
be at equilibrium, and there may be species where this alternative historical interpretation of the 
model is plausible. 

We have also modeled migration between populations as occurring at single, instantaneous time 
points. This is, of course, a dramatic simplification of the migration process. This model will work 
best when gene flow between populations is restricted to a relatively short time period. Situations 
of continuous migration violate this assumption and lead to unclear results (Supplementary Figure 
14). The relevance of this assumption will depend on the species and the populations considered. 
In humans, the relevance of continuous versus discrete mixture events is an open question-some 
aspects of genetic variation appear compatible with continuous migration [61], while other aspects 
do not [37]. Indeed, both sorts of models are likely relevant at different time scales [62]. 

We have also rely on the implicit assumption that the history of the species being analyzed is 
largely tree-like. We have made this assumption to simplify the search for the maximum likelihood 
graph; additionally, we speculate that in graphs with complex structure, there will be many graphs 
that lead to identical covariance matrices, and thus several different histories will be compatible 
with the data. That said, improvements to the search algorithm could allow the assumption of 
approximate treeness to be somewhat relaxed. Currently, if the number of admixed populations is 
large relative to the number of unadmixed populations, this assumption breaks down. For example, 
in the human data, note that we see no evidence of the documented gene flow from Neandertals 
to all non- African populations [39] (Figure 3B). The reason for this is that the large number 
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of populations with admixture can be accommodated in the tree by allowing the branch from 
Neandertals to Africans to be slightly underestimated (additionally, by using SNPs ascertained in 
Africa, we have selected against sites that are informative about Neandertal ancestry). If only a 
single non- African population is included in the analysis, the relationship between Neandertals and 
the non- African population is clearer (Supplementary Figure 15). 

Conclusions. A number of extensions to the sort of model described here are of potential in- 
terest. First, the historical relationships between populations could be useful as null demographic 
models for the detection of natural selection [48, 63, 64]. Second, in a given individual, the best-fit 
graph relating the individual to other populations may change along a chromosome; this sort of 
information could be of use in local ancestry inference. Finally, we have not used the information 
about demographic history present in linkage disequilibrium; approaches that explicitly use this 
information may provide additional power to detect migration events and estimate their timing, at 
an additional computational cost [20, 53, 65]. 

Methods 

Graph estimation. As described in the Results, we developed an algorithm called TreeMix 
that uses the composite likelihood in Equation 28 to search for the maximum likelihood graph. 
Estimation involves two major steps. First, for a given graph topology, we need to find the maximum 
likelihood branch lengths and migration weights. Second, we need to search the space of possible 
graphs. First consider a given graph topology. We iterate between optimizing the branch lengths 
and weights. If the edge weights are known, the observed entries of the covariance matrix can be 
written down as an overdetermined system of linear equations (as in Equations 13-15). We solve 
this system by non- negative least squares [66]. Though the least squares solution is the maximum 
composite likelihood solution in the case where all entries of the covariance matrix have equal 
variance, it is not strictly the maximum likelihood solution in cases with unequal variances. The 
algorithm could be extended to unequal variances using a weighted least squares approach, but we 
have not implemented this. We then do a golden section search for the optimal weight (between 
zero and one) on each migration edge [67]. At each step in the golden section search, we update the 
branch lengths. We optimize the weight of each migration edge in turn, and iterate over migration 
edges until convergence. 

To search the space of possible graphs, we take a hill-climbing approach. We start by finding 
a local optimum tree, taking an algorithmic approach similar to Felsenstein [30]. We randomly 
select three populations, optimize the branch lengths for all three possible trees, and choose the 
best (in terms of the composite likelihood) tree. Then, we add the remaining populations one by 
one in a random order. To add a population, we try attaching it to all branches of the current 
tree, optimizing the branch lengths for each one as described above, and find the most likely spot. 
We then perform a round of local rearrangements (i.e., nearest-neighbor interchanges [50]) around 
each internal node, keeping the resulting tree only if it increases the likelihood. 
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After adding all populations, we calculate the residual covariance matrix, R. We then add 
migration edges in a directed matter. First, we find the M pairs of populations with the maximum 
residuals (these are the pairs of populations with the worst fit under the model). In the results 
reported, M = 4. We define a "neighborhood" around each population of a pair as the tips within 
a distance of E edges of the focal population. In applications above, we use E = 3. This defines a 
set of pairs of populations that either have a poor fit, or are located in the graph near populations 
with a poor fit. We take each of these pairs in turn. For each pair, we identify the set of nodes 
in the path from each member of the pair to the root of the graph. This gives us two sets of 
nodes. We take all pairwise combinations of nodes in each set, and look at residuals between the 
populations that are the descendants of each node. If all of the residuals are positive, we add a 
migration edge between the two nodes and estimate its maximum likelihood weight. We then keep 
only the single edge that most increases the likelihood of the graph. After adding a migration edge, 
we attempt nearest-neighbor interchanges at the source and destination of the migration event, 
attempt changing the source and destination of all migration events, and attempt changing the 
direction of all migration arrows. Once we have reached the local maximum by this method, we 
attempt nearest-neighbor interchanges at all internal nodes. We iterate over this procedure for a 
predetermined number of migration edges. We then test the migration edges for significance as 
described. 

The TreeMix source code is available at http://treemix.googlecode.com. 

Three- and four-population tests of treeness. We implemented three- and four-population 
tests as described in Reich et al. [37]. For the relationship between the /—statistics and the 
covariance model underlying TreeMix, see the Supplementary Material. For the three-population 
test, we estimated as in Reich et al. [37], and tested whether is it less than zero. We report the 
Z-score for this test. To obtain a standard error on the estimate of f^, we used a block jackknife 
similar to Reich et al. [37]. However, Reich et al. [37] split the genome into blocks based on 
distance (with variable numbers of SNPs per block); we split the genome into blocks of K SNPs 
(and thus the blocks will be of variable size) . 

For the four-population test for treeness, we calculate the statistic as in Reich et al. [37], 
and test whether it is different than zero. Again, we report a Z-score for this test. Standard errors 
for the /4 statistic were obtained as for the fa statistic. 

Human data. The human data we used were downloaded from http : //www . cephb . f r/ en/hgdp/ 
on August 16th, 2011 (the data set labeled Harvard HGDP-CEPH genotypes). They consist of 
several panels of SNPs ascertained from low-coverage genome sequencing of single individual from 
different populations and then genotyped in the Human Genome Diversity Panel [54]. Additionally, 
at each site, a single sequencing read from the Denisova and Neandertal genome sequencing projects 
was sampled and the allele reported. These data have the property that they allow for complete 
control of the ascertainment strategy, and allow us to test the robustness of inference to different 
ascertainment schemes. For the main analyses, we used the panel of autosomal SNPs ascertained 
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in a single Yoruban individual; there are 124,115 such sites. For some analyses, we also used the 
panel of autosomal SNPs ascertained in a single French individual; there are 111,970 such sites. For 
all analyses with TreeMix, we used a window size (-K) of 500; this corresponds to a window size 
of approximately 10Mb. For all TreeMix analyses, we set the Neandertal and Denisova samples as 
the outgroups. 

Since we have only a single allele from the Neandertal and Denisova populations, we cannot 
calculate heterozygosity in these populations for unbiased estimation of the covariance matrix 
(see Supplementary Information). To account for this, we simply chose a relatively low level of 
heterozygosity and assigned it to both populations. In the Yoruba ascertained SNPs, we used 
a heterozygosity of 0.13, and for the French ascertained SNPs, we used a heterozygosity of 0.2. 
In practice, this only affected the lengths of the terminal branches to Neandertal and Denisova; 
running TreeMix with a heterozygosity of zero in both populations resulted in the same graph 
topologies (not shown). 

Dog data. Allele counts for the dog breeds and wild canids reported in Boyko et al. [58] were 
downloaded from http://genome-mirror.bscb.cornell.edu/ on July 30, 2011. These data con- 
sist of counts of reference and alternate alleles at 61,468 sites in 85 dog breeds and wild canids. We 
removed the Jackal and Scottish Deerhound for having relatively high amounts of missing data, 
and the village dogs because it is unclear if they represent a coherent population. We also removed 
all SNPs on the X chromosome. This left us with 60,615 SNPs in 82 populations. We ran TreeMix 
with a window size (— K) of 500. This corresponds to a window size of approximately 20 Mb. For 
all TreeMix analyses, we set the coyote as the out group. 

The ascertainment scheme used for SNP discovery in dogs was complicated [68] . The largest set 
of SNPs were ascertained by virtue of being different between the boxer and poodle assemblies. This 
should lead to an overestimation of the distance between the boxer and the poodle in our analysis. 
Indeed, in Figure 5B, a considerable negative residual between the boxer and poodle is visible. 
Another set of SNPs were ascertained by being heterozygous within a boxer individual, and a third 
set were ascertained by comparison between a boxer and wild canids. These latter SNPs should 
lead to an overestimation of the distance between the boxer and the wolf in our analysis (as we see 
for the poodle); in fact, we infer migration between the boxer and the wolf. This ascertainment 
issue may have led us to underestimate the amount of gene flow in the comparison. 

Simulations. All simulations were performed using ms [69]. The exact commands used are listed 
in the Supplementary Material. When running TreeMix on simulations without ascertainment, we 
used a window size of 5000 SNPs; for simulations with ascertainment we used windows of 1000 
SNPs. Consensus trees were generated using SumTrees v. 3. 1.0. [70] 
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Correcting covariances for finite sample size. In the main text, we define the variance- 
covariance matrix W of allele frequencies between populations without accounting for sampling 
variance. Here, we show the calculations corrected for sample size. Consider n biallelic loci typed 
in m populations of diploid individuals, and let the sample size in population i at locus k be 
(with missing data, the number of individuals can vary across loci). Let the counts of the two alleles 
in population i at locus k be rii k and 2Ni k — rii k (with one allele being arbitrarily defined as the 
reference in all that follows), the true allele frequency in the population be Xik, and the observed 
allele frequency be = j^r- We assume the rii k are binomially distributed with parameters 
2N ik and X ik , and are independent for all i and k. Recall that the allele frequency in the ancestral 
population is xa, and that the covariance between populations i and j with respect to the ancestral 
frequency xa is V^. We begin by defining Vy using the observed allele frequencies at a single SNP 
k: 

\ lj = E[{X lk -x A ){X 3k -x A )] (1) 
= E[[(X lk - X ik ) + (X ik ~ xa)] [{X jk - X jk ) + {X jk - x A )]\ (2) 
= E[(X lk - x A )(X jk - x A )} + E[(X lk - X ik ){X jk - X jk )}. (3) 

The bias in the estimate of is thus E[(X ik — X ik ) 2 } if i = j (i.e., it is the sampling variance in 
X ik ) and zero otherwise. This follows from the fact that the rii k are assumed to be independent 
across i. 

Now consider all n SNPs, and let the mean bias across all SNPs be B{. At a given SNP k, 
the sampling variance in population j is X ik is Xtk ^N-f' tk ' > (from the binomial sampling of Xi k ), 
so the mean bias across SNPs is proportional to X ik (l — X ik ) (i.e., the mean across all SNPs of 
Xi k (l — Xik)). A natural estimator of Bi is then: 

B - = w, < 4 » 

where hi is an unbiased estimate of the heterozygosity in population i averaged over all SNPs [Nei, 
1978]: 

1 ^ n ik (2Nj - n ik ) 
hi ~n^ N t (2N t - 1) • (5) 

As derived in the main text, the sample covariance of populations i and j, Wjj, is: 

1 m ^ m ^ m m 

w v = - ^ E v ^ - - E v ^ + -2 E E v ^'- 

k=l k=l k=l k'=l 

The bias in the estimate of Wy (let us call this B[-) is then: 

Tjl T R Bi Bj Y^k = l B k r „x 

rfii - I [i=j]Bi 1 5 \ 7 > 

j mm m z 
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where I[i=j] is an indicator that evaluates to 1 if i = j and zero otherwise. We can then estimate 
the unbiased covariance Wa as: 





(8) 



n 



where Hk = = ^ lh ■ If there is missing data in either population i or population j, we simply 
ignore the SNP for that pairwise comparison of populations. Since the mean allele frequency across 
populations is important here, large amounts of missing data (or correlated missingness between 
populations) could result in skewed covariances. We thus exclude populations with large amounts 
of missing data. 

Nonidentifiability of the drift parameters in an admixed population. In the main text, 
we write down a model for the allele frequencies in an admixed population, and claim that the 
amount of genetic drift occurring before and after the mixture event are nonidentifiable. Consider 
the graph in Supplementary Figure 1. We can write down the expected variances and covariances 
involving the admixed population: 



and we are interested in estimating w, c±, C2, and C3. It is clear from the above that c±, C2, and 
C3 do not appear except as a linear combination. Adding additional populations does not add 
additional information about these parameters, unless they are assumed to result from the same 
mixture event. 

We choose to set C2 and c\ to zero, and estimate only C3, which can now be thought of as a 
composite branch length that sums all the three components of genetic drift. A subtle point is that 
all of this drift is weighted by (1 — w). When estimating w, then, the true relative contributions of 
ci, C2, and C3 could lead to a bias in the estimation of w. For example, if c\ and/or C2 are large, this 
could bias the estimation of (1 — w) upwards. We believe this is likely the cause of the downward 
bias in w in the simulations in Figure 2D in the main text. 

Graph representation of the TreeMix model. In the main text, we describe a specification 
of V (the variance/covariance matrix of allele frequencies, defined with respect to an ancestral 
population) in terms of a system of linear equations. A useful alternate notation describes V in 
terms of a graph [Roller and Friedman, 2009]. Let G be a rooted, directed, acyclic graph with a set 
of nodes N and a set of directed edges E. Each edge e has an associated length, c e , and a weight, 
w e (between zero and one). A special class of edges, called migration edges, are forced to have 
length zero. The sum of weights of edges entering a given node is one. There is one node which is 
the root (a node with only outgoing edges), and each population corresponds to a tip (a node with 



V12 



V 22 
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w c 5 x A [l - x A ] 

[ci + w 2 (c 2 + c 5 ) + (1 - w) 2 (c 3 + c4)]z A [l - x A ] 
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only incoming edges). 

Define {P{} to be the set of all possible paths in G from the root to the tip corresponding to 
population i (if the graph is a tree, there is only one such path). Each individual path p has a 
weight, w{p) = H e £ P w e . Now define the overlap between two paths as: 

0(pi,Pj) = w(pi)w(pj)I[e e pj]c e (12) 

where I[e € Pj] is a function that evaluates to one if edge e is in pj, and zero otherwise. We can 
now write down the expected covariance between populations i and j as: 

V ^'= E E ( 13 ) 

Pi€{Pi}pje{Pj} 

In the special case where G is a tree, there is only one path per population and all of the edges have 
weight one, and so Vjj reduces to a sum of the lengths of branches shared by the two populations. 

Relationship of this model to /— statistics. Tests for "treeness" in three and four-population 
trees [Keinan et al., 2007; Reich et al., 2009] have used a framework in which the distances between 
populations are quantified in terms of "/—statistics" comparing the allele frequencies between 
the populations. Below, we briefly describe these tests in the notation of our model. Consider 
the expected fa statistic calculated between populations 1, 2, and 3, with corresponding allele 
frequencies Xi, X2, and X 3 . 

fa(X 1 ;X 2 ,X 3 ) = E[(X 1 -X 2 )(X 1 -X 3 )} (14) 
= E[[(X! - x A ) - (X 2 - x A )] [{Xi - x A ) - (X 3 - x A )}] (15) 
= Vn- V12-V13 + V23. (16) 

Consider the situation where populations 1 and 3 form a clade relative to 2 (i.e., population 2 is 
an outgroup). If population X\ is not admixed, this reduces to: 

fa(X 1 ;X 2 ,X 3 ) = V 11 -V 13 . (17) 

This is necessarily greater than zero (since V13 <= Vn). If X\ is admixed, then V12 can be 
important and the fa statistic can be negative. A test for a negative fa statistic is thus a test for 
admixture in population X\ [Reich et al., 2009]. However, this signal can be weakened by large 
amounts of drift in X\ (i.e., a large Vn), or mixture between X 2 and X 3 [Reich et al., 2009]. 
Similarly, consider the expected fa statistic computed on the tree [[1,2], [3, 4]], where 1, 2, 3, and 
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4 are populations, and Xi, X2, X% and X4 are the corresponding allele frequencies: 



f i (X 1 ,X 2 ;X 3 , X 4 ) = - X 2 )(X 3 - X*)] 

= - x A ) - (X 2 - s A )][(*3 - x A ) - (Xi - x A )}] 



(18) 
(19) 
(20) 
(21) 



= V13 - V 23 - V14 + V 24 



If the tree is correct (i.e., if populations 1 and 2 are a clade relative to populations 3 and 4), all of 
these quantities are zero. A test for a non-zero fi statistic is thus a test for treeness [Reich et al., 



Simulation commands. For all simulations, we used ms [Hudson, 2002]. To generate the tree- 
like data depicted in Figure 2A in the main text, the command is: 

ms 400 400 -t 200 -r 200 500000 -I 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 

20 20 20 20 20 -en 0.00270 20 0.025 -ej 0.00275 20 19 -en 0.00545 19 0.025 -ej 0.00550 

19 18 -en 0.00820 18 0.025 -ej 0.00825 18 17 -en 0.01095 17 0.025 -ej 0.011 17 16 -en 
0.01370 16 0.025 -ej 0.01375 16 15 -en 0.01645 15 0.025 -ej 0.01650 15 14 -en 0.01920 
14 0.025 -ej 0.01925 14 13 -en 0.02195 13 0.025 -ej 0.02200 13 12 -en 0.02470 12 0.025 
-ej 0.02475 12 11 -en 0.02745 11 0.025 -ej 0.02750 11 10 -en 0.03020 10 0.025 -ej 0.03025 
10 9 -en 0.03295 9 0.025 -ej 0.03300 9 8 -en 0.03570 8 0.025 -ej 0.03575 8 7 -en 0.03845 
7 0.025 -ej 0.03850 7 6 -en 0.04120 6 0.025 -ej 0.04125 6 5 -en 0.04395 5 0.025 -ej 
0.04400 5 4 -en 0.04670 4 0.025 -ej 0.04675 4 3 -en 0.04945 3 0.025 -ej 0.04950 3 2 

-en 0.05220 2 0.025 -ej 0.05225 2 1 

To create trees with considerably longer branch lengths, we multiplied all branch lengths by 50: 
ms 400 400 -t 200 -r 200 500000 -I 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 

20 20 20 20 20 -en 0.135 20 0.025 -ej 0.1375 20 19 -en 0.2725 19 0.025 -ej 0.275 19 
18 -en 0.41 18 0.025 -ej 0.4125 18 17 -en 0.5475 17 0.025 -ej 0.55 17 16 -en 0.685 
16 0.025 -ej 0.6875 16 15 -en 0.8225 15 0.025 -ej 0.825 15 14 -en 0.96 14 0.025 -ej 
0.9625 14 13 -en 1.0975 13 0.025 -ej 1.1 13 12 -en 1.235 12 0.025 -ej 1.2375 12 11 
-en 1.3725 11 0.025 -ej 1.375 11 10 -en 1.51 10 0.025 -ej 1.5125 10 9 -en 1.6475 9 
0.025 -ej 1.65 9 8 -en 1.785 8 0.025 -ej 1.7875 8 7 -en 1.9225 7 0.025 -ej 1.925 7 
6 -en 2.06 6 0.025 -ej 2.0625 6 5 -en 2.1975 5 0.025 -ej 2.2 5 4 -en 2.335 4 0.025 
-ej 2.3375 4 3 -en 2.4725 3 0.025 -ej 2.475 3 2 -en 2.61 2 0.025 -ej 2.6125 2 1 

For simulations with migration, we added a migration event approximately 100 generations be- 
fore the present. For example, migration from population 1 to population 10: 

ms 400 400 -t 200 -r 200 500000 -I 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 
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20 20 20 20 20 -em 0.002675 10 1 4000 -en 0.00270 20 0.025 -em 0.00270 10 1 -ej 0.00275 
20 19 -en 0.00545 19 0.025 -ej 0.00550 19 18 -en 0.00820 18 0.025 -ej 0.00825 18 17 
-en 0.01095 17 0.025 -ej 0.011 17 16 -en 0.01370 16 0.025 -ej 0.01375 16 15 -en 0.01645 
15 0.025 -ej 0.01650 15 14 -en 0.01920 14 0.025 -ej 0.01925 14 13 -en 0.02195 13 0.025 
-ej 0.02200 13 12 -en 0.02470 12 0.025 -ej 0.02475 12 11 -en 0.02745 11 0.025 -ej 0.02750 
11 10 -en 0.03020 10 0.025 -ej 0.03025 10 9 -en 0.03295 9 0.025 -ej 0.03300 9 8 -en 
0.03570 8 0.025 -ej 0.03575 8 7 -en 0.03845 7 0.025 -ej 0.03850 7 6 -en 0.04120 6 0.025 
-ej 0.04125 6 5 -en 0.04395 5 0.025 -ej 0.04400 5 4 -en 0.04670 4 0.025 -ej 0.04675 
4 3 -en 0.04945 3 0.025 -ej 0.04950 3 2 -en 0.05220 2 0.025 -ej 0.05225 2 1 

For simulations of populations exchanging migrants on a lattice, we used the following command: 

ms 200 1500 -t 40 -r 40 100000 -I 10 20 20 20 20 20 20 20 20 20 20 -ma x 
x 0.1 0.1 0.1 0.1 x 0.1 0.1 0.1 0.1 0.1 x 0.1 0.1 
0.1 0.1 x 0.1 0.1 0.1 0.1 0.1 0.1 0.1 x 0.1 0.1 0.1 0.1 0.1 0.1 
0.1 x 0.1 0.1 0.1 0.1 x 0.1 0.1 0.1 0.1 0.1 x 0.1 
0.1 0.1 0.1 x -ej 0.41 10 1 -ej 0.41 9 1 -ej 0.41 8 1 -ej 0.41 7 1 -ej 0.41 6 1 -ej 
0.41 5 1 -ej 0.41 4 1 -ej 0.41 3 1 -ej 0.41 2 1 



Discussion of simulation errors. In Figure 2C in the main text, we showed that TreeMix was 
extremely accurate in most simulation situations. However, there are a few situations in which it 
performed poorly. Most notably, this was for simulated admixture between population 1 and 5. The 
errors in these simulations tended to be of the same type (Supplementary Figure 7). Additionally, 
in the simulations of migration from population 15 to population 20 with a weight of 10%, there 
was also a considerable error rate. However, these errors were not consistent across simulations, 
and are likely due to the algorithm simply not detecting the admixture event at all. 

Analysis of human data including Oceanian populations. As described in the main text, in 
the human HGDP data we used two sets of allele frequencies with different ascertainment schemes- 
one at SNPs ascertained by sequencing a single Yoruban individual, and one at SNPs ascertained by 
sequencing a single French individual. We initially ran TreeMix on both data sets using all popula- 
tions to estimate the maximum likelihood trees. The trees estimated using the two ascertainment 
schemes are nearly identical (Supplementary Figure 11). We then used TreeMix to identify migra- 
tion events. The algorithm arrived at quite different conclusions about the Oceanian populations 
in the two different data sets (recall that these are the exact same individuals, just genotyped 
at different SNPs) (Supplementary Figure 15). In the Yoruba-ascertained data, the East Asian 
populations are inferred to be admixed, with the Melanesians as a source population. However, 
in the French-ascertained data, the Oceanians are inferred to be admixed. When the Oceanian 
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populations are excluded from analysis, the algorithm comes up with nearly the same graph in 
both datasets (Figure 4 in the main text and Supplementary Figure 10). 

It is not immediately clear why there is a discrepancy between these two datasets when looking at 
Oceanian populations. However, Oceania has a particularly complicated genetic makeup, involving 
at least four distinct components of ancestry: Denisovan gene flow, Neandertal gene flow, native 
Oceanian, and gene flow from Austronesian speakers [Reich et al., 2010, 2011; Wollstein et al., 
2010]. These different components of ancestry may be picked up to differing extents by SNPs from 
the different ascertainment panels, leading to conflicting results. 

List of migration events inferred in the human data. Here we list the ten migration edges 
inferred in the human data and present in Figure 4 in the main text: 

1. Yoruba ->■ Mozabite, w = 19% 

2. Mandenka — > (Palestinian, (Bedouin, Mozabite)), w = 13% 

3. Mandenka — > Druze, w = 6% 

4. Mankenka — > (Brahui, Makrani), w = 6% 

5. Ancestral non- African — > Cambodian, w = 17% 

6. (All Europe, Middle East) Hazara, w = 46% 

7. (All Europe, Middle East) Uygur, w = 46% 

8. All Native Americans — > Russian, w = 12% 

9. (Tuscan(French(Italian, (Basque, Sardinian) — > Maya, w = 12% 
10. Mozabite — > (Tuscan(French(Italian, (Basque, Sardinian), w = 22% 

List of migration events inferred in the dog data. Here we list the ten migration edges 
inferred in the dog data and present in Figure 6 in the main text: 

1. (Greyhound, Whippet) — > Borzoi, w = 47% 

2. (AlaskanMalamute,SiberianHusky) — > Samoyed, w = 47% 

3. Wolf Basenji, w = 25% 

4. (ChowChow,ChineseSharPei) ->■ (Pekingese, Shih Tzu), w = 28% 

5. Bulldog Bull mastiff, w = 31% 

6. Boxer — > Chinese shar-pei, w = 8% 

7. Siberian Husky — > (Akita,(ChowChow,ChineseSharPei)) w = 17% 
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Figure 1: A graph with a mixture event. Capital letters represent nodes, branch length 
parameters are in blue, and weight parameters are in red. 

8. Wolf -»■ Boxer, w = 8% 

9. Mastiff— > Saint Bernard, w = 13% 

10. Whippet — > Italian Greyhound, w = 37% 
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Figure 2: Replicates of inferred trees from simulated data. We generated tree-like data 
using the topology in Figure 2 A in the main text. In figure 2B in the main text, we show the 
inferred tree with mean branch lengths. In this figure, we show four representative individual trees. 
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Figure 3: Inferred trees on ascertained data. We generated tree-like data using the topology 
in Figure 2 A in the main text. We then used only the SNPs that were polymorphic in either 
population 1 (A.) or population 20 (B.) to infer the trees. The correct topology was obtained in 
all 100 simulations; the branch lengths in each figure are the mean across all simulations. 
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Figure 4: Histogram of p-values for migration in simulated data. We generated 100 tree- 
like datasets using the topology in Figure 2A in the main text. We then randomly chose two 
populations (without replacement), added a migration edge between the two populations, and 
tested for significance using the procedure described in the main text. Plotted is the histogram of 
p-values for the significance test. If the p-values are properly calibrated, this distribution should be 
uniform (dotted line). Though the distribution is not completely uniform, there is no skew towards 
low p-values. 
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A. Consensus tree 
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Figure 5: Consensus tree in simulations with long branches. We generated 100 tree-like 
datasets using the topology in Figure 2 A in the main text, multiplying all branch lengths by 50. 
We then inferred the maximum likelihood tree. A. Plotted are the mean branch lengths from 
the simulations. All simulations resulted in the same inferred topology. B. In each simulation, 
we scaled the residuals by the average standard error, then averaged these scaled residuals across 
simulations. Plotted are the mean scaled residuals across the 100 simulations. The most extreme 
residuals are not large (around 0.3 standard errors), but tend to be present between closely related 
populations. 
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A. Example tree 1 
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B. Residual fit from tree in A. 



pop1 




T-WCO^^tDNffiO)OT-OICO^IOCDN(DO)0 
Q.CLCLQ.Q.CLQ.Q.CLT-T--i-t-T--i-T-T--i--i-(M 

o o o o o o o o o aaaaaaaaaaa 

Q.Q.Q.Q.Q.Q.Q.Q.CLOOOOOOOOOOO 

aaaaaaaaaaa 



C. Example tree 2 
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D. Residual fit from tree in C. 
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Figure 6: Example trees from simulations with long branches. We generated 100 tree-like 
datasets using the topology in Figure 2 A in the main text, multiplying all branch lengths by 50. 
We then inferred the maximum likelihood tree. In Supplementary Figure 5, we show the average 
inferred tree. Here, we show two representative trees (A. and C. and the residuals corresponding 
to each tree (B. and D. 
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A. correct graph topology 
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B. common error (w = 1 0%) 
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Figure 7: Representative errors in simulations. We examined the simulations in which 
TreeMix did not reach the correct answer. A. The correct topology for the simulations presented 
in the other panels. B. A representative example of an incorrect topology inferred from the sim- 
ulations of a migration event with weight 10% from population 1 to population 5 (this topology 
accounted for all observed errors). C. A representative example of an incorrect topology inferred 
from the simulations of a migration event with weight 30% from population 1 to population 5 (this 
topology accounted for 95% of all errors). 
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Figure 8: Replicate graphs inferred in the human data. These graphs were generated in 
an identical manner as Figure 4 in the main text, but using different random input orders for 
populations during tree-building. All random input orders gave very similar results. 
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Figure 9: Residual fit from graph of human data presented in the main text. Plotted are 
the residuals from the fit of the graph presented in Figure 4 in the main text. 
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Figure 10: Graph inferred from SNPs ascertained in a single French individual. The 

graph was generated in an identical manner as Figure 4 in the main text, but using a panel of 
SNPs ascertained in a single French, rather than a single Yoruban, individual. The inferred graph 
is extremely similar to that in Figure 4. The one major difference is that, in this graph, the Mozabite 
appear as an admixture of a Sardinian population rather than a Middle Eastern population; this 
configuration is seen in some runs of TreeMix on the Yoruba-ascertained data (Supplementary 
Figure 8A) 
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A. ML human tree (Yoruba) 



B. ML human tree (French) 
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Figure 11: Trees inferred using the human data including the Oceanians. We show the 
maximum likelihood trees and residuals for the human data including the Oceanian populations, 
plotted in the same manner as in Figure 3 in the main text. Trees were inferred using the panel of 
SNPs ascertained in a single Yoruban individual (A. and C.) and the panel of SNPs ascertained 
in a single French individual (B. and D.). See Supplementary Material for discussion. 
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A. Yoruba ascertainment 
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Figure 12: Graphs inferred using the human data including the Oceanians. We show the 
maximum likelihood graphs for the human data including the Oceanian populations, plotted in the 
same manner as in Figure 4 in the main text. Six migration edges were inferred in each graph. 
Graphs were inferred using the panel of SNPs ascertained in a single Yoruban individual (A.) and 
the panel of SNPs ascertained in a single French individual (B.). See Supplementary Material for 
discussion. 
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Figure 13: Residual fit from graph of dog data presented in the main text. Plotted are 
the residuals from the fit of the graph presented in Figure 6 in the main text. 
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A. Simulated grid 
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Figure 14: TreeMix run on populations with continuous migration. We simulated a set of 
populations on a lattice, where each population has constant gene flow at a rate of ANm = 0.1 with 
neighboring populations. All populations split from an outgroup 16,400 generations in the past. 
The exact ms command is given in the Supplementary Material. The configuration of the lattice 
is presented in A. (population 1 is the outgroup). TreeMix inferred no consistent tree structure. 
Three representative trees are presented in B.-D. 
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A. Africa tree B. Africa tree residuals 
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Figure 15: TreeMix run on human data using only a single non-African population. 

We inferred the maximum likelihood tree (A.) using only the African populations and one non- 
African population (French), using SNPs identified in a single Yoruban individual. In examining 
the residuals (B.), a relationship between the French and the Neandertal is clear. We then inferred 
three migration events (C), where we do see that the French contain some Neandertal ancestry 
(w = 1.2%). Residual fit for this graph is shown in D. 
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